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Abstract 

We consider the solution of complex symmetric shifted linear sys- 
tems. Such systems arise in large-scale electronic structure simulations 
and there is a strong need for the fast solution of the systems. With 
the aim of solving the systems efficiently, we consider a special case of 
the QMR method for non-Hermitian shifted linear systems and propose 
its weighted quasi-minimal residual approach. A numerical algorithm, 
referred to as shifted QMR_SYM(_B), is given by the choice of a particu- 
larly cost-effective weight. Numerical examples are presented to show the 
performance of the shifted QMR_SYM(B) method. 

1 Introduction 

In this paper we consider the solution of complex symmetric shifted linear sys- 
tems of the form 

(A + a e I)x w =b, £=l,2,...,m, (1) 

where A(ai) :— A + ogl are nonsingular N-by-N complex symmetric sparse 
matrices, i.e., A(at) = A(ai) T ^ A(ai) T , with scalar shifts ae <E C, / is the N- 
by-N identity matrix, and b are complex vectors of length N. The above 
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systems arise in large-scale electronic structure simulations [T5], and there is a 
strong need for the fast solution of the systems. 

Since the given shifted linear systems |T]) are a set of sparse linear systems, 
it is natural to use Krylov subspace methods, and moreover since the coefficient 
matrices are complex symmetric, one of the simplest ways to solve the shifted 
linear systems is applying (preconditioned) Krylov subspace methods for solving 
complex symmetric linear systems such as the COCG method [16], the COCR 
method [13] , and the QMR_SYM method [3] to all of the shifted linear systems 
(fT]). On the other hand, denoting the rt-dimensional Krylov subspace with 
respect to A and b as K n (A, b) := span{6, Ab, . . . , A n ~ 1 b}, we observe that 

K n (A,b) = K n (A{at),b). (2) 

This implies that once basis vectors are generated for one of the Krylov sub- 
spaces K n (A(ae), b), these basis vectors can be used to solve all the shifted lin- 
ear systems. In other words, there is no need to generate all Krylov subspaces 
K n {A(at), b), and thus computational costs involving the basis generation, e.g., 
matrix- vector multiplications, are saved. Here we give a concrete example: if we 
apply the conjugate orthogonal conjugate gradient (COCG) method to all the 
linear systems (JTJ, then bases for K n (A(ag),b) are generated for I = 1, 2, . . . , m. 
On the other hand, if we apply the COCG method to just one of the shifted 
linear systems (J]) (referred to as the " seed system"), then the Krylov basis 
vectors are generated from the seed system and these vectors are used to solve 
the rest of the shifted linear systems. 

Based on the observation ([2]), the shifted COCG method [15] has been re- 
cently proposed for solving complex symmetric shifted linear systems. The 
feature of the shifted COCG method is that the method performs COCG on a 
seed system and makes it possible to complete COCG for all shifted linear sys- 
tems without further matrix- vector multiplications. The feature is completely 
different from some of the other well-known shifted linear solvers such as the 
shifted BiCGStab(£) method [6J or the restarted shifted GMRES method [7J 
since these perform BiCGStab(^) (or GMRES) on a seed system but a different 
method on the shifted linear systems in order to keep the residuals colinear. 
The feature of the shifted COCG method plays a very important role in the 
seed switching technique [14] that avoids a minor problem of the shifted COCG 
method: one requires the choice of a seed system and an unsuitable choice may 
lead to the undesirable result that some shifted linear systems remain unsolved. 

There is another approach to solving the shifted linear systems (JTJ. That 
is the use of Krylov subspace methods for non-Hermitian shifted linear systems 
such as the shifted BiCGStab(£) method [6], the shifted (TF)QMR method 
[4 , the restarted shifted FOM method [10] , and the restarted shifted GMRES 
method [7], see also, e.g., [11]. We readily see that the relation ([2]) holds not 
only for complex symmetric matrices but also for non-Hermitian matrices, and 
these methods are based on the use of this shift-invariant relation. Therefore, 
this can be a good approach. However, since these methods do not exploit the 
property of complex symmetric matrices, their computational costs can be more 
expensive than that of the shifted COCG method. 
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In this paper we consider the shifted QMR_SYM method that is a special 
case of the QMR method for non-Hermitian shifted linear systems [4] and clar- 
ify the most time consuming part of it for a large number of shifted linear sys- 
tems. Then, in order to reduce the cost, we propose a weighted quasi-minimal 
residual (WQMR) approach and propose a specific weight. We experimentally 
show the practical efficiency of the resulting algorithm, referred to as shifted 
QMR_SYM(B), when the number of shifted linear systems is large enough. 

The present paper is organized as follows: in the next section, shifted 
QMR_SYM is described in order to specify the most time consuming part for a 
large number of shifted linear systems. In Section 3, we propose a WQMR ap- 
proach with a specific weight for reducing the cost of the most time consuming 
part. The resulting algorithm, shifted QMR_SYM(£?), and its properties are 
given. In Section 4, some results of numerical examples from electronic struc- 
ture simulations are shown to see the performance of the shifted QMR_SYM(i3) 
method. Finally, some concluding remarks are made in Section 5. 

Throughout this paper, unless otherwise stated, all vectors and matrices are 
assumed to be complex. M, M T , M H = M denote the complex conjugate, 
transpose, and Hermitian transpose matrix of the matrix M. \\v\\w denotes the 
W-norm written as (v^Wv) 1 / 2 , where W is Hermitian positive definite. 

2 The QMR.SYM method for solving complex 
symmetric shifted linear systems 

In this section, the shifted QMR_SYM method and its properties for solving 
complex symmetric shifted linear systems are introduced. 

The QMR method for shifted linear systems was first formulated in [4] for 
the case of a general non-Hermitian matrix. Therefore, by simplifying the non- 
Hermitian Lanczos process [5] , as is known from other papers such as [31 116j , a 
shifted simplified QMR method, shifted QMR_SYM, is readily obtained for the 
case of a complex symmetric matrix. Its algorithm is given next. 

Algorithm 2.1. Shifted QMR.SYM 

«S° = p { -\ = p ( o } = o,-i - V(b T b) 1/2 , g[ e) = (b T b) 1/2 , 

for n = 1, 2, . . . do: 

(The complex symmetric Lanczos process) 
a n = v^Av n , 

V n +1 = AV„ - a n V n - fln-lVn-l, 
Pn = (vl +1 V n+1 ) 1/2 , 
V n +1 = Vn+l/Pn, 

t {e) - 3 , tW - a +<j e t {e) - 3 

(Solve least squares problems by Givens rotations) 

for I = 1,2, ... , m do: 
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if ||rW|| 2 /||6|| a >e ) then 
for i = max{l, n — 2}, . 



»+l,n , 




M 

(£) _ z n+l,n (i) 
°n M) °n i 

t w - 



(Update approximate solutions cc« ) 

f'n \ L n-2,nl L n-2,n-2)* J n -2 \ t 'n-l,n/ v n-l,n-l)t'n-l ' 

™W - a-W + r W /*(*) } D W 
end if 
end 

if ||r^|| 2 /||6||2 < e for all I, then exit. 

end 

Algorithm 2.1 can be regarded as a natural combination of the results given in 

HE]- 

In order to know that the numerical solution is accurate enough, one may 
need to compute the residual 2-norms. In that case, the following computation 
may be useful. 

Proposition 1 (See |5j) The 2-norms of the nth residuals of the approximate 
solutions Xn of the shifted QMR_SYM method are given by 

Wrtfh = l&Hk3-illa for £ = l,2,...,m, 

i {£) (i) (I) . U) , (I) 

where — — s„ 'Wn + c„ v n+ i and w\' = V\. 

Proposition 2.1 is a result known to hold for the QMR method [5]. Therefore, 
it also holds for the above specialized variant. The rest of this section describes 
some special properties of the shifted QMFLSYM method. 

Proposition 2 (See |2j) Let A € R JVxAr be real symmetric, ag € C be com- 
plex shifts, and b 6 R N . Then the shifted QMRJSYM method (Algorithm 2.1) 
enjoys the following properties: 
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/. all matrix-vector multiplications can be done in real arithmetic; 

H. an approximate solution at nth iteration step for each I has minimal 
residual 2-norms, i.e., x$ 's are generated such that min||rn over 
x$ eK n (A,b); 

®- \\r$h = \&\ for 1 = 1, 2,..., m, n > 0. 

The above properties are known results since the properties have been proved 
for each individual shift. See [5] for details. The properties of Proposition 2.2 
may be very useful for large-scale electronic structure simulations [15j and a 
projection approach for eigenvalue problems [12) since there are complex sym- 
metric shifted linear systems to be solved efficiently under the assumption of 
Proposition 2.2. 

3 An iterative method for solving complex sym- 
metric shifted linear systems 

In this section we consider complex symmetric shifted linear systems with a 
large number of shifts. For such systems, say m 3> 1, the most time consuming 
part of Algorithm 2.1 can be generating approximate solutions since the cost 
for the recurrences is 6Nm + 3m per iteration step. In order to reduce the cost, 
in this section we propose a weighted quasi-minimal residual approach with a 
specific weight. In the next subsection we discuss the details of the approach. 
In Subsection 3.2 we give a specific weight to achieve the reduction of the cost. 

3.1 A weighted quasi-minimal residual approach 

Before the formulation of a Weighted Quasi-Minimal Residual (WQMR) ap- 
proach, let us recall the complex symmetric Lanczos process (see, e.g., Algo- 
rithm 2.1 in [3]). 

Algorithm 3.1. The complex symmetric Lanczos process 

set fa = 0, v = 0, r ^0e C N , 
set vx = r /(rjr ) 1/2 , 
for n = 1, 2, . . . , m — 1 do: 
a n = v„Av n , 

v n +i = Av n - a n v n - /3 n _it> n _i, 

&n = (vl +1 Vn+l ) 1/2 , 
V n +1 = V n+ i/(3 n . 

end 

The matrix form of Algorithm 3.1 is known as follows: let T, i+ i,„ and T n be 
the (n + 1) x n and n x n tridiagonal matrices whose entries are recurrence 
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coefficients of the complex symmetric Lanczos process, which are given by 
fax [3! \ 
Pi ct 2 



V 



'•• Pn-1 

Pn-l OL n 

Pn I 



T n := 



fax Px 

Pi a 2 



\ 



\ 



'■■ Pn-l 
Pn-1 a n I 



and let V n be the N x n matrix with the Lanczos vectors as columns, i.e. 
V n := (vx, V2, ■ ■ ■ , v n ). Then from Algorithm 3.1, it follows that 



AV n = V n+x T n+ x,n = V n T n + p n v n+1 el, 



(3) 



where e„ = (0,0, ...,1) T G R". 

Now we are ready to describe the WQMR approach. Let x„ ' be approximate 
solutions at the nth iteration step for the systems Jl}, which are given by 



x 



w - v v w e - i 2 m 



(4) 



where y„ £ C™. Then, from the definition of residual vectors Tn '■= b — 
(A + <j&I)x n , , the update formulas ((4|), and the matrix form of the complex 
symmetric Lanczos process ([3]), we readily obtain 



(5) 



r$ = K+i (giex - Ti% <n yW) , where T^ 1>n := T n+1 , n + *t 

Here, e x is the first unit vector written by ei=(l, 0, . . . , 0) T and gx = (6 T 6) 1 ^ 2 . 
It is natural to determine yiP such that all residual 2-norms ||ri^||2 are mini- 

(£) 

mized. However, this choice for y n is impractical due to large computational 
costs. Hence, an alternative approach is given, i.e., the vectors y n are deter- 
mined by solving the following weighted least squares problems: 



yj^' = are min 
z<«ec« 



(6) 



where W n +x is an (n + l)-by-(n+ 1) nonsingular matrix. Thus W^ +1 W n +x can 
be used as a weight since it is Hermitian positive definite. One of the simplest 
choices for W n +x is the identity matrix. In this case, from W„+i = In+i we 
have 



arg mm 



yi^i J n , ln ^ n 



(7) 



The vector that is minimized is called quasi-residual. Algorithm 2.1 is obtained 
by solving using Givens rotations, see, e.g., p. 215]. 

A slightly generalized choice proposed in [5] is W n +i = '■— diag(a>i, u>2, ■ ■ 

w n +i) w hh cjj > for all i. Then, it follows that 



y 



arg mm 



uxgiex - n n+1 T^l l n z^ 
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Of various possible choices for uii, a natural one is u>i = \\vi\\2 since fi n +i 
contains then the diagonal entries of the upper triangular matrix R n +i that is 
obtained by the QR factorization of V n +i- If we choose W n +i = R n +i, where 
Vn+i — Qn+iRn+i, then from §5§ and ^ we have 



mm 

z%>ec» 



5i e i 



r W At) 

n+l,n n 



= mm 
= min 
= min 



giRn+iei - Rn+iT„ +lin z$ 

Qn+iR n +i(giei -Tnl l>n z^ 

V n+1 (g iei -T^ ltn zW 
r W\\ 

n II 2' 



By solving the above weighted least squares problems, all residual 2-norms are 
minimized. Hence W n +i = £l n +i is a rational choice. For each individual shift, 
the resulting algorithm is the same as that in [3l Algorithm 3.2]. 

3.2 A choice of the weight suitable for a large number of 
shifts 

In the previous subsection, we have described the WQMR approach and men- 
tioned that the choice of the weight W^jM^+i with W n +i = I n +i leads to 
the shifted QMRJ3YM method (Algorithm 2.1). Under the assumption of 
proposition 2.2, the shifted QMR_SYM method is ideal in the sense of Faber- 
Manteuffel's theorem [T] since it enjoys minimal residual property and requires 
not long-term but short-term recurrences for updating approximate solutions, 
and thus one may think that there is no need to choose other possible weights. 
However, we will show in this subsection that even under the assumption of 
proposition 2.2 there is a practical weight for the WQMR approach. The moti- 
vation for the choice of the weight mainly comes from the freedom of the number 
m of complex symmetric shifted linear systems. 

Now we consider the computational costs of Algorithm 2.1 for a large number 
of complex symmetric shifted linear systems, i.e., m 3> 1. For the case m 3> 1 
we readily see from Algorithm 2.1 that computing the recurrences for updating 
approximate solutions is the most time-consuming part due to a cost of 6Nm + 
3m per iteration step. Hence we will now consider a weight to reduce the 
computational cost for the recurrences of x„ . To achieve this we propose the 
following choice: 

W n+1 = L% such that = ( of ) . (8) 
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where B n is an n-hy-n upper bidiagonal matrix of the form 



/At) M 
I c i,i 1 



B (t) 



At) 

L 2.2 



V 



At) 

L n-l.n 
t W J 



(i) 

and is lower triangular and will be specified below. 

Next we derive recurrence formulas for updating the approximate solutions 
Xn\ From © with the choice W n+ \ = of © it follows that 



,{t) 



arg mm 
arg min 



arg mm 

2< f) ec« 



(9) 



yl-^n+l^l ^ri+l J n+l,n z n 



B (0 



2 \ y n +i 



r (<) 



From the above least squares problems we readily see that y„ = (-Bn ) 1 g^f' > - 

Hence it follows from (HJ) and using (p 1 p 2 ■ ■ ■ p n ) := ^(-B^-*) -1 that we have 
the following coupled two-term recurrence relations: 



X 



(v -t {e) v {e) )/t w 

n ' = X n—1 9n ^ Pn ^ ■ 



The cost per iteration is now 5Nm. Substituting p\ = t\-p\ into the above 
recurrences, we have the even more efficient recurrence formulas. 

Pn ^ u n k s n— l,n/ c n-l,n-l )Pn-li 



n-1 + (ffri V^n,n)jPn ^ • 



By this reformulation, the cost becomes 4iVm + 2m. 

Next, let us consider a choice for satisfying (JSJ). Let i 7 ^ be an n-by-n 

matrix for the form 



(I, 



\ 



if 1 



(10) 



and let T be tridiagonal. Then, by determining such that the + i) entry 
of the matrix F n l) (i)T is zero, we can fulfill ||SJ| in the following way: 



At) 



At) 



B 



(t) 



(11) 
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From the above we see that ^+i(n)-Fi+i(n - 1) • • • = L^{ 1; and thus 

and L^n are related by 

4'|i=^iW(o? °) forn = 2,3,..., (12) 

where = ^2^(1)- F rom the above we readily see that the matrices are 
lower triangular with all l's on the diagonals, and this property will be used for 
the proof of proposition 3.1 given later. The shifted QMFLSYM method with 
the weight (L^ ] +1 ) H L { r p +1 is referred to as shifted QMR_SYM(S). Its algorithm 
is given next. 



Algorithm 3.2. Shifted QMR.SYM(B) 



4 ] = P { -\ = P ( o ] =0,v 1 = &/(& T &) 1/a > gf> = (b T 6)V 2 , 
for n = 1,2,... do: 

(The complex symmetric Lanczos process) 

a n = v„Av n , 

v n+ i = Av n - a n v n - (3 n -iv n -i, 

Pn = {vl +1 V n+ i) 1/2 , 
V n +1 = V n+1 //3 n , 

t n ~l,n — Pn-li tn}n = a n + &ti ^n+l,n ~ @ni 

(Solve weighted least squares problems) 
for I = 1,2, ... , m do: 
if ||rW|| 2 /||6|| 2 >£, then 

for i = max{l, n— l},...,n — 1 do: 

M _ MM , M) 

L i+l,n ~ 1% L i,n ' L i+l,n> 

end 

M) 

M) _ b n+l.n 

Jn It) i 

t (e) - 

„W - M-xW 

i/n+l — Jn y-n ) 

(Update approximate solutions xffi) 

— \ L n-l,nl i n-l,n-l)l'n-li 

X W - £C W + Co W /t^ )T) (e) 
^n ~ A n-1 ' Kiln I l n,n)l'n i 

end if 
end 

if ||rW|| 2 /||6|| 2 < e for all £, then exit. 
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end 



Similar to Proposition 2.1, there is an efficient way to evaluate residual 2-norms 
as follows: 

(£) 

Proposition 3 The nth residual 2-norms of the approximate solutions x n for 
the shifted QMR_SYM(B) method are given by 

\\r$h= \^U-\\v n+1 h for e=l,2,...,m. 

Proof. The proof is similar to that of Proposition 2.1. It follows from ([5]), 
©, ©, and recalling = (B^)' 1 ^ that we have 

From (fTOf and (jlip is a lower triangular matrix with all l's on the di- 

agonals. Thus, 1 is also a lower triangular matrix with all l's on the 

diagonals. It follows that (L^^) 1 e n+ i = e n+ i, and thus we have r«' = 

3n+i^n+ie„+i = .g^j^n+i, which concludes the proof. 

When we solve a large number of shifted systems, the computational cost 
of the residual 2-norms for the shifted QMR_SYM(_B) method is much cheaper 
than the that for the shifted QMR_SYM method since the former cost is of order 
N and the latter is of order Nm per iteration step. 

Observing algorithms of the shifted QMR.SYM and the shifted QMR_SYM(B) 
methods, we see that the work for the weighted least squares problems and 
for updating approximate solutions in the shifted QMR_SYM(i3) method is 
lower than that in the shifted QMR_SYM method. For the case m ^> 1, 
the most time-consuming part is updating the approximate solutions. In this 
case, the shifted QMR_SYM(i?) method may be more efficient than the shifted 
QMR_SYM method since the shifted QMR_SYM(B) method requires 4Nm+2m 
operations per iteration step for the update while the shifted QMR_SYM method 
needs 6Nm + 3m operations. This difference will be clearer when we use the 
results of Propositions 2.1 and 3.1 as stopping criterion. On the other hand, in 
terms of number of iterations, under a certain assumption, the convergence of 
the shifted QMR_SYM(B) method is worse than that of the shifted QMR.SYM 
method but not worse than that of the shifted COCG method, as described by 
the following proposition: 

Proposition 4 Under the assumption of Proposition 2.2, the shifted QMR_SY 
M(B) method (Algorithm 3.2) enjoys the following properties: 

I. all matrix-vector multiplications can be done in real arithmetic; 

H. if breakdown does not occur and each matrix T n -\-o~iI n is nonsingular, then 
^ r (j),SQ(B)^_^ r (e),scocG^^ r (e),SQ^ g _ i2,...m, where the su- 
perscripts SQ(B), SCOCG, and SQ are short for shifted QMR_SYM(B), 
shifted COCG, and shifted QMRSYM, respectively; 
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m. \\r^' SQ(B) \\ 2 = \g^l 1 \for£=l > 2,...,m, n>0. 

Proof. The proof of (I) is the same as that of Proposition 2.2, and is based 
on the fact that under the assumption the complex symmetric Lanczos process 
generates real basis vectors. 

Next, we give a proof of (II). The nth residuals of the shifted COCG 
method [TS] belong to b — (A + a£l)K n (A + agl,b). Hence, each r W' SCOCG 

, ., , (£),SCOCG ■ / A . r\ t r («,SCOCG , t r • , -i 

can be written as r„ = o — (A + ael)V n yn , where V n is the 

same matrix as in ((3|) Since each r ^' SCOCG [ s orthogonal to each subspace 
K n (A+a e I,b){=K n (A,b)}, i.e., r (f)> SCOCG jl K n (A,b), wehaveV^b-V^(A+ 
aeI)V n y W > SCOGG = 0, and thus it follows from © that we have the rela- 
tion y W< SCOCG = {V J(A + a e I)V n }^b = g 1 (T^e 1 , where := 
T n + a e I n . Since the shifted QMR_SYM(B) method has the form ©, it is suf- 
ficient to show that y W> SCOCG = y^^ B \ From © and JT2JI it follows that 
y^ SQ ^ = (B?)- 1 ^ = ff i(^)- 1 [4|0]4'i 1 [er|0] T = 9i(B^)-^ ei = 
Since from (|lip and (|12|) we can readily confirm the re- 
lation Ln^Tn^ — Bn \ we have yn? ,S ™ B '* — gi(Tn^)~ 1 ei, which is the same as 
j/(^),SCOCG. The inequality in (77) follows from Proposition 2.2 since under the 
given assumption the shifted QMFLSYM method enjoys the minimal residual 
property. 

Finally, we give a proof of (HI) . If follows from the proof of (I) that \\vi ||a = 1 
for all i. Thus from Proposition 3.1 we have \\rn ^' SO -( B ^ || 2 = |<?„i_i| • ||vn+i||2 = 

Ifl&l for*=l,2,...,m, n>0. 

In the property (II) of Proposition 3.2, breakdown may occur due to the 
choice (|10[) of the weighted least squares problems. 

From proposition 3.2 we see that in terms of the number of iteration steps the 
shifted QMR_SYM(B) method never converges faster than the shifted QMR_SYM 
method but it converges at the same iteration step as the shifted COCG method 
does. Since the efficiency of the shifted COCG method has been shown already 
and for the case m 3> 1 the computational cost of the shifted QMR_SYM(i?) 
method is much less than that of the shifted QMR_SYM method, the shifted 
QMR_SYM(S) method can also be useful. This will be supported by some 
numerical examples in the next section. 



4 Numerical examples 

In this section, we report on some numerical examples with the shifted COCG 
method, the shifted QMR_SYM method (Algorithm 2.1), and the shifted QMR_SYM(B) 
method (Algorithm 3.2). We evaluate these methods in terms of computa- 
tional time. All tests were performed on a workstation with a 2.6GHz AMD 
Opteron(tm) processor 252 using double precision arithmetic. Codes were writ- 
ten in Fortran 77 and compiled with g77 -03. In all cases the stopping criterion 
was set as e = 10 -12 . 
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3.1. Example 1 



The first problem comes from the electronic structure calculation of a bulk 
Si(OOl) with 512 atoms in [15] . which is written as follows: 

(ail-H)xW =ei, £=l,2,...,m, 

where a t = 0.400 + (£-l + i) /1000, H E R 2 °^x 2 °^ is a symmetric matrix with 
139264 entries, e x = (1,0, . . . ,0) T , and m = 1001. Since the shifted COCG 
method requires the choice of a seed system, we have chosen the optimal seed 
(£ = 714) in this problem; otherwise some linear systems will remain unsolved. 

Figure 4.1 shows the number of iterations of each method to solve the £th 
shifted linear systems. For example, in Fig. 4.1, the number of iterations for 
the shifted COCG method at I = 600 is 150, which means the shifted COCG 
method required 150 iterations to obtain the (approximate) solutions of the 
600th shifted linear system, i.e., {a^I — iJ)a;( 600 ' = e\. 




200 400 600 800 

Serial number of shifted linear systems (I) 



Figure 1: Number of iterations for the shifted COCG method, the shifted 
QMR_SYM method, and the shifted QMFLSYM(B) method versus the serial 
number of the shifted linear systems. 

From Fig. 4.1 we make three observations: first, the three methods required 
almost the same number of iterations at each l\ second, in terms of number of 
iterations, the shifted QMFLSYM method often converged slightly faster than 
the other two methods. This phenomenon is closely related to Proposition 2.2, 
which will become clearer later; third, for each method the required number 
of iterations depends highly on the shift parameters an. This result may come 
from varying eigenvalues of the coefficient matrices ail — H since if we choose 
ai close to an eigenvalue of H, then ogl — H is close to singular. Conversely, 
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Figure 2: Logio of the relative residual 2-norms versus the number of iterations 
of the shifted COCG method, the shifted QMFLSYM method, and the shifted 
QMR_SYM(B) method for the shifted linear system with t = 701, i.e., 0701 = 
1.100 + 0.00K 

from the shape of Fig. 4.1 one may obtain a partial distribution of eigenvalues 
of if. 

One of the residual 2-norm histories for the three methods is given in Fig. 4.2. 
From Fig. 4.2 we see that the relative residual 2-norm of the shifted QMR_SYM 
method decreases monotonically and at every iteration step the norm is less 
than those of the other two methods. Hence we can say that the property 
(if) of Proposition 2.2 is experimentally supported by this history. We make 
another observation in the histories of the shifted COCG method and the shifted 
QMR_SYM(B) method: During the first about fifty iterations, the two methods 
show the same histories. After that their histories varies gradually. Hence we 
see that the property (if) of Proposition 3.2 is also experimentally supported 
by these histories. 

Each computational time of the three methods is given in Fig. 4.3, where the 
horizontal axis denotes the number m of shifted linear systems that are solved. 
For example, in Fig. 4.3, the computational time of the shifted COCG method 
at m = 200 is about 0.76 [sec], which means that it required about 0.76 [sec] 
to solve the shifted linear systems: {(0.400 + 0.001i)i - H}x^ = ei, {(0.401 + 
O.OOli)/ - H}x^ = ei, . . . , {(0.599 + 0.001i)7 - H}x {2m ) = ei . From Fig. 4.3 
we see that as the number m grows larger, the shifted QMFLSYM method 
required more CPU time than the other two methods. On the other hand, 
the shifted QMR_SYM(f?) method required almost the same CPU time as the 
shifted COCG method. This phenomena can be attributed to the computational 
costs of updating approximate solutions for each method since there are three 
facts: first, we know from Fig. 4.1 that the three methods required almost 
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Figure 3: Required CPU time given in seconds versus the number of shifted 
linear systems for each iterative method. 
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Figure 4: The ratio of each computational time to the one of the shifted COCG 
method versus the number of shifted linear systems. 
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the same iterations; second, the shifted QMR_SYM(i?) method requires almost 
the same computational costs as the shifted COCG method, while the shifted 
QMRJSYM method tends to require more computational costs per iteration 
than the other two methods; third, for large m, updating approximate solutions 
is one of the most time-consuming parts. The second and the third facts have 
been discussed already in the previous two sections. 

In Fig. 4.3 we can see little about the properties of the three methods for 
small I. We therefore show the ratio of each computational time to the computa- 
tional time of the shifted COCG method in Fig. 4.4. We see from Fig. 4.4 that in 
terms of the ratio of CPU times the shifted QMR.SYM method and the shifted 
QMR_SYM(_B) method converged much faster than the shifted COCG method 
when the number of shifted linear systems is small, say, m < 200. This can be 
explained in the following way: for small m, updating approximate solutions 
does not affect the CPU time so much. Other operations such as matrix-vector 
multiplications are now the most time-consuming parts since the three methods 
required almost the same number of iterations, see Fig. 4.1. From Proposition 
2.2 (7) and Proposition 3.2 (I) we know that in this case the shifted QMRJSYM 
method and the shifted QMR_SYM(_B) method require real matrix-real vector 
multiplications. On the other hand, the shifted COCG method requires real 
matrix-complex vector multiplications. Moreover, dot products and vector ad- 
ditions of the complex symmetric Lanczos process used in the shifted QMR_SYM 
method and the shifted QMR_SYM(_B) method can be done in real arithmetic. 
Hence, the two methods converged much faster than the shifted COCG method. 

3.2. Example 2 

The second problem comes from the electronic structure calculation of bulk 
fee Cu with 1568 atoms in 45] and is given as follows: 

(ad - H)x^ = ei, £ = l,2,...,m, 

where a t = -0.5 + (I - 1 + i)/1000, H E #14112x14112 is a symmetric ma t r i x 
with 3924704 entries, e x = (1, 0, . . . , 0) T , and m = 1501. 

Each computational time of the three methods for solving the m shifted 
linear systems is shown in Fig. 4.5, where horizontal and vertical axes are the 
same as those used in Fig 4.3. The ratio of each computational time to the 
computational time of the shifted COCG method is shown in Fig. 4.6. From 
Figs. 4.5, 4.6 we see that the size of this matrix is about 7 times larger than 
before, and the three methods show similar tendencies to those seen in the 
previous example. 

5 Concluding remarks 

In this paper, the shifted QMR_SYM method was described as a specialization 
of the QMR method for general non-Hermitian shifted linear systems [4]. The 
method has an advantage over the shifted COCG method in that there is no 
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Figure 5: Required CPU time given in seconds versus the number of shifted 
linear systems for each iterative method. 
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Figure 6: The ratio of each computational time to the one of the shifted COCG 
method versus the number of shifted linear systems. 
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need to choose a suitable seed system. On the other hand, we have found that 
for a large number of shifted linear systems, the most time-consuming part of 
the shifted QMFLSYM method is updating approximate solutions, and this cost 
is higher than that of the shifted COCG method. We therefore have proposed 
the weighted quasi-minimal residual approach with a specific weight for reduc- 
ing the computational cost for updating approximate solutions. The resulting 
method, shifted QMR_SYM(S), also does not require to choose a suitable seed 
system, which is an advantage over the shifted COCG method. From numer- 
ical experiments we have learned that shifted QMR_SYM and QMR_SYM(S) 
are competitive in comparison to the shifted COCG method. In particular, 
QMR_SYM(B) can be the method of choice for solving complex symmetric 
shifted linear systems with a large number of shifts that arise from large-scale 
electronic structure theory. In future work, numerical tests for general complex 
symmetric shifted linear systems will be done to evaluate the performance of 
the method. 
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